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ABSTRACT  distinguish  it  from  source  and  channel  vocoders, 

and  to  suggest  that  identification  of  formants  and 
In  tills  paper  we  make  the  point  that  a  wide  other  mode  parameters  proceeds  statistically, 

variety  of  spectrum  types  admit  to  modal  analysis  The  domain  of  attraction  for  our  modal  decon- 

wherein  the  modes  are  characterized  by  amplitudes ,  position  of  the  covariance  sequence  includes  ARMA 

frequencies,  and  damping  factors.  The  associated  processes,  harmonically-  or  nonharmonically-related 

modal  decomposition  is  appropriate  for  both  con-  sinusoids,  white  noise,  and  linear  combinations  cf 

linuous  and  discrete  components  of  the  spectrum.  these.  This  means  the  decomposition  generalizes 

The  domain  of  attraction  for  the  decomposition  the  covariance  sequence  model  implicit  in  DFT, 

includes  ARMA  sequences,  harmonically-  or  nonhar-  Pisarenko  [3],  and  linear  prediction  formulations 

monical Ly-related  sinusoids,  damped  sinusoids,  of  the  parametric  spectrum  analysis  problem  in  the 

white  noise,  and  linear  combinations  of  these.  same  way  that  Prony's  ancient  model  [1]  general- 

Numcrlcal  results  are  presented  to  illustrate  izes  DFT,  almost  periodic,  and  autoregressive 

the  identification  of  mode  parameters  and  corro-  models  for  data, 

sponding  spectra  from  finite  records  of  perfect 

and  estimated  covariance  sequences.  The  results  COVARIANCE  SEQUENCE  APPR0X1MANT 

for  sinusoids  and  sinusoids  in  white  noise  are 


interpreted  in  terms  of  inphase  and  quadrature 
effects  atLri Durable  to  the  finite  record  length. 


INTRODUCTION 


Motivated  by  the  form  of  the  covariance 
sequence  for  line  and  rational  spectra,  we 
propose  the  following  order-p  covariance  sequence 
approximant : 


This  is  a  paper  about  parametric  covariance  P  t 

sequence  approximation  for  the  purposes  of  model  rt(p)  “  ^  Aizi  *  t*0,l,2,... 

identification  and  spectrum  analysis.  As  with  all  i“l 


parametric  methods  of  analysis  the  key  idea  is  to 

pose  a  defensible  model  for  the  data,  or  some  r_t(p)  *  rt(P) 

important  statistical  descriptor  such  as  the  covar¬ 
iance  sequence,  and  then  to  identify  model  pa mm- 

eters.  z ^  “  p^e 

Our  approach  Is  to  represent  the  covariance 


sequence  of  a  wlde-sense  stationary  process  in  a 
modal  decomposition.  Amplitudes,  frequencies 
(formants),  and  damping  constants  then  become  mode 
parameters  to  be  identified.  These  mode  parameters 
may  be  used  in  vibration  analysis,  speech  analysis, 
and  forecasting  and  control. 

Corresponding  to  the  modaL  decomposition  for 
the  covariance  sequence  is  a  modal  decomposition 
for  the  underlying  process.  Thus  the  decomposition 
becomes  un  anal j sis/synthesis  tool  applicable  to 
data-compressed  speech  and  data  communication,  and 
to  stochastic  simulation.  The  resulting  synthe¬ 
sizer  might  he  termed  a  statistical  vocoder  to 

’supported  by  the  Office  of  Naval  Research  under 
Contract  N00014-75-C-0518  and  the  Army  Research 


We  call  7}.  the  i  1  complex  mode,  with  frequency  - . 
and  radius  p^,  and  A^  the  corresponding  node  weignt. 
The  weights  and  modes  appear  in  complex  conjugate 
pairs.  Typically  a  subset  of  the  p^  are  unity  to 
model  the  discrete  component  of  the1  sped rum,  and 
the  rest  are  less  than  unity  to  model  the  continu¬ 
ous  part  of  the  spectrum.  See  [21  for  review  of 
discrete-time  spectral  theory. 

The  Prony  Device:  Begin  with  a  finite  covari¬ 
ance  string  (rg,r^f . . . ,r ^ , . . . , )  ,  assumed  to 
be  the  first  2p  Fourier  Scries  coefficients  for  e 
density  f(u>).  Solve  for  the  regression  coeffi¬ 
cients  (n^t<«2«'*ata  )  that  satisfy 
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P-1 


P 

rp+l 


P-2 

r2p-2  rp  Vl  ||“p|  L  ‘2p-i: 

Form  the  poJynonial 

P(z)  *  i  +  a,z  *  +  ...  +  a  z  p 
1  P 

and  find  its  roots  zj*z2»'*“»z  .  Construct  the 
Vandermonde  patrix  P 


..P-  i 
*1 


P"1 


and  solve  for  weights  A  ,A  , . 
ing  system  of  equations: 


. , A  in  the  follov- 
P 


A, 
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1 
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r  . 

Pw 

P-1J 

The  modal  decomposition 
P  t 

r  (p)  '  Z  A  z 
i-1 


reconstructs  (r  )  and  extends  it  as  a 

linear  conhimtl Ion  of  complex  modes.  If 
r2  j)»  P  >  Pjj,  comes  from  an  ARMA  (Pg»Pg)  prucess 
then  A  ^  "  0  for  p  >  and  the  modal  decomposition 
matches  r ^  for  all  t.  Here  <  p^. 

The  Fcjer  Device:  The  trigonometric  sum 


f(w) 


l  r  (p)p 1 
t—p 


j  tw 


ANALYSIS  AND  SYNTHESIS 

Assume  we  are  given  a  modal  decomposition  of 
the  form 


rt(p) 


S  Az" 
1-1  1 


The  following  question  naturally  arises:  is  there 
a  finite-dimensional  spectrum  representation  (or 
vocoder)  that  generates  stochastic  realizations  of 
the  wide-sense  stationary  process  x£  that  has  rc(p) 
for  its  covariance?  The  answer  Is  a  qualified  yes, 
and  the  development  goes  as  follows. 

Let  X  be  a  complex  process  generated  by  the 
autoregressive  scheme 

xt  -  pVi  +  £t 

2 

where 

XQ  :  N(0,A) 

Et  :  sequence  of  i.i.d.  N(0 , (1- | p | A) 
random  variables 


This  process  has  covariance 
t 


Rt  -  A  p‘ 


t-0,1,... 


The  process 


Jut 


xt  ■  V 

with  u  fixed  and  independent  of  X£  has  covariance 


Az 


pe 


ju 


Think  of  x  generated  this  way  as  x  (A,p,u).  Gen¬ 
erate  p  independent  processes  like  this  of  the  fcnv 
xt  -  x(.(Ai  ,p1>ui)  and  sum  to  get 


i-1 

with  covariance 


P 

Z  A, 


i-1 


Ju. 


This  statistical  vocoder  is  illustrated  In  Figure  1. 
It  generates  a  special  class  of  processes  with 
modal  decompositions  such  that  A^  >  0.  It  captures 
DFT  and  Pisarenko  decompositions  as  special  cases. 


uniformly  approximates  f(u)  as  ptl  and  pv».  Note 
the  term  simply  draws  the  modes  z|C^  into 

II  1 

(pZjT  -corresponding  to  pole  locations  more  Inte¬ 
rior  to  the  unit  disc. 


The  notation  X„:N(0,A)  Indicates  X  is  complex 
normal  with  mean  zero  and  varlancc°A>0.  The  model 
la  trivially  generalized  to  other  distributions. 


PARAMETER  IDENTIFICATION  FROM  A  FINITE  RECORD 


modified  squared  error 


In  this  section  we  propose  a  two-step  (a-one, 
a-two)  procedure  for  obtaining  a  modified  least 
squares  fit  of  r  (p)  to  a  covariance  string  that 
has  been  estimated  as  follows; 


r 

t 


N 

r 

i«l 


ii+t 


E 

m 


N-l 

Z 

t-p 


P 

(  Z 
1-0 


ao 
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This  amounts  to  filtering  the  errors  with  a  moving 
average  filter  and  counting  them  after  p  steps 
iiave  elapsed.  If  the  filter  weights  are  selected 
so  that 


Least  Squares:  Define  the  squared  error 
between  rt(p)  and  r^  as  follows: 

N-l 

E  -  I  (r  -r  (p))z 

t— '.'+1 


Z  a  z 
i=0 


-i 


It  (1-2.2  l) 
1-1 


then  the  approximating  r  (p)  sequence  satisfies 
this  homogeneous  difference  equation  on  its  tail: 


Setting  derivatives  of  E  with  respect  to  (A^,z  ), 
i*l, 2, to  zero  yields'  a  discrete  form  of  the 
Aigrain-Wil 1 equations  [6],  slightly  modified 
to  reflect  the  fact  we  are  approximating  a  two- 
sided  sequence: 


JE_ 

3A, 


il-1 

Z 

t~  M+l 


=  0  (1-1,2,...  ,p) 


tlx  airt-i(p) 


t  >  p 


This  is  one  of  the  things  Prony  saw  while  studying 
the  effects  of  alcohol  vapor  pressures  in  1795. 

With  this  modification  the  contribution  of 
r^(p)  to  the  modified  squared  error  is  annihilated 
and  we  are  left  with 


JLL 

3z . 


iJ-1 

Z 

t*-D+l 


(rt-rt(p))Ai|t|zi 


t  -1. 


d-1.2 . P) 


Using  the  syixr.etry  of  the  and  r^fp)  sequences 
we  may  write  this  system  of  equations  as 

PTQr  -  PTqPA 

PTQr  »  P1QPA 


where 


■r 


(ru  V 


Zp-H 
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Q- 


2  0 
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0 

0  2 


"rN-l> 


lA . V 

P  j 

I  A,  P 
i-1  1  3Zl 


This  system  of  equations  is  nonlinear  in  the  param¬ 
eters  (A  ,z.)  and  calls  for  iterative  procedures 
if  the  solution  is  to  be  "exact".  We  suggest  here 
a  modification  of  these  equations  for  which  a 
tractable  solution  procedure  exists. 

A  Modification  f»  Solution:  Consider  the 


N-l  P  7 

E  -  Z  (  Z  a  r  Y 
m  t-p  i-0  1  C  1 

Minimisation  with  respect  to  the  a.  leads  to  a 
covariance-method  of  linear  prediction  on  the  tail 
of  the  covariance  sequence, 

T  T 

R  Ra  -  -R  r  , 


where 


P-1 


N-l  N-2 


-r-  -  (V Vl . rN-l> 


M-P 


~P 

T 


(*!•••. *ap) 


The  elements  of  R  R  arc  "covariances  of  covari¬ 
ances". 

There  are  a  variety  of  fast  algorithms  for 
solving  these  equations.  Once  the  are  found 
the  corresponding  modes  are  obtained  using  a 
polynomial  rootfinding  routine  such  as  Muller's 
algorithm  [A].  With  these  modes  determined  we  may 
solve  a  linear  system  of  equations  for  the  mode 
weights  Aj.  This  is  basically  Prony's  method  at 
work  except  for  the  fact  that  the  ai  are  obtained 
from  a  squared  error  criterion  rather  than  Prony’o 
point-wise  criterion  that  results  f  om  sotting 
N-2p: 


4 


Ra  =  r 

One  cannot  help  wondering  what  path  statisti¬ 
cal  inference  might  have  taken  had  Gauss  been 
aware  of  Prony's  parametric  interpolation  formu¬ 
lae  when  he  published  his  work  on  least-squares 
in  1809.  For  at  this  early  date  all  the  ingred¬ 
ients  of  a  least-squares  theory  of  rational 
approximation  would  have  existed. 

APPLICATIONS  AND  NUMERICAL  RESULTS 

Here  we  present  the  results  of  several  numeri¬ 
cal  experiments.  These  experiments  are  summar¬ 
ized  in  Table  1.  In  the  table  the  column  labels 
nave  the  following  meanings: 

Exact  Model-Underlying  model  that  generated 
data 

No.  Samples  -  Number  of  samples  used  to  esti¬ 
mate  r^  (-»  means  exact  covariance  used) 

No.  Corr.  Lags  -  Number  of  correlation  lags 
used  in  fitting  algorithm 

Fitted  ARMA( •  ,  • )  -  For  the  modal  decomposi¬ 
tion  the  initial  order  denotes  the  initial 
number  of  poles  identified.  When  the  final 
order  differs  from  the  initial  order  this 
means  a  noise  pole  (at  z*0)  was  added.  For 
the  modified  least  squares  (MLS)  the  Initial 
order  tells  what  order  AR  was  fit  to  ini¬ 
tialize.  The  final  order  gives  the  approx¬ 
imating  AKMA . 


A  100 

rt  “  Too  1^1  xi*i+t 

See  Figure  4  for  spectra  of  approximating  modal 
decomposition  and  approximating  modified  least 
squares  fit. 

two  Closely-Spaced  Sines:  Here 

x  ■  sin  7-  t  +  cos. 95  7-  t 
t  4  4 

160 

rt  °  160  1f1  xixi+t 

Sec  Figure  5  for  spectra  of  approximating  modal 
decomposition  and  approximating  modified  least 
squares  fit. 

ARMA(1 ,1) :  Here  x  is  the  output  of  the 
following  filter  excited  by  white  noise: 

h(2)  .  idL!5£i 

1-0. 25z 

See  Figure  6  for  spectra  of  approximating  modal 
decomposition  and  approximating  modified  least 
squares  fit. 

ARMA(3,2) :  Here  x  is  the  output  of  the 
following  filter  excited  by  white  noise: 


Two  Closcly-S|)ac.ed  Sines  In  WCN:  Here 

r  •  cos  7  t  +  0.01  cos. 95  ;  t  +  100  6 
c  4  4  t 


H(z) 


l-1.75z~1+0.8z  2 _ 

l-1.5z_1+1.21z'2-0.455z'3 


See  Figure  2  for  spectrum  uf  approximating  modal 
dccomposlt ion . 

Sine  in  AH  Noise:  Here 

*  cos  it  t  ♦  Z  *{H(z)H(z  *)  } 

H(z)  -  l/(l-J0.8z"1)(l-t-j0.8z'1) 

See  Figure  3  for  spectrum  of  approximating  modal 
decomposition. 

Comment :  In  a  1976  paper  on  digital  filter  design 

1 6 J  one  of  us  advocated  a  two  step  procedure  for 
fitting  long  AR  sequences  to  exact  correlation  se¬ 
quences,  followed  by  modified  least  squares  fitting 
of  an  AKMA(p,q).  The  modified  least  squares  pro¬ 
cedure  was  originally  proposed  by  Kalman  and  subse¬ 
quently  fully  developed  by  Mullis  and  Roberts  [7]. 
At  the  end  of  that  paper  we  suggested  that  the 
same  procedure  might  be  applied  to  AKMA  spectrum 
analysis.  This  suggestion  we  have  explored  in  eon- 
} unction  with  our  modal  th*compos 1 1 ion  studies. 
Spectral  res u Lis  for  modified  least  squares  filling 
from  long  AR  models  are  shown  Logethcr  with  speelra 
for  modal  decompositions  in  the  following  examples. 

Sine:  Here 


See  Figure  7  for  spectra  of  approximating  modal 
decomposition  and  approximating  modified  least 
squares  fit. 

More  examples  may  be  found  in  [5]. 

CONCLUSIONS 

We  have  proposed  an  approach  to  covariance 
sequences  and  rational  spectrum  approximation  that 
captures  discrete  spectra  and  ARMA  spectra  as  spe- 
c'al  cases.  The  approach  is  based  on  a  modal  de¬ 
composition  for  the  covariance  sequence  wherein 
the  modes  correspond  to  poles.  For  a  special  clas 
of  processes  with  modal  decomposition  there  is  a 
stochastic  synthesis  algorithm  (or  statistical  vo¬ 
coder)  that  may  be  used  to  generate  realizations. 
This  could  be  a  valuable  algorithm  for  data-com- 
pressed  communication  based  on  the  analysis/syn¬ 
thesis  ideas  of  this  paper. 

Numerical  results  for  noise-free  and  noisy 
sinusoids  are  encouraging.  A  virtue  of  this 
technique-and  a  virtue  we  want  to  emphasize-is 
that  one  obtains  from  the  analysis  technique  both 
n  spectrum  and  a  map  of  the  underlying  modes  of 
the  covariance  sequence.  These  modes  may  he  mon¬ 
valuable  in  some  instances  than  the  spectrum  it¬ 
self.  In  fact,  the  spectrum  f(u>)  can  often  ob¬ 
scure  very  interesting  fine  structure  in  the  dat.* 
that  one  can  observe  directly  in  the  covariance 
sequence  npproximant  r  (p). 
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Igure  2.  Spectral  estimate  for  2  sines  in  white 
noise.  Given  15  known  covariance 
lags.  SNR  =  -20  dB,  -40  dB 


Figure  1.  Statistical  Vocoder 
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Figure  4.  Spectral  catiratc  for  sine.  Using  20  bias- 
estimatej  la:  s  from  100  samples 
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Figure  5.  Spectral  estimate  for  two  closely-spaced 

sines.  Using  20  bias-estimated  covariance 
lags  from  160  samples 
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Figure  6.  Spectral  estimates  for  ARMA(1,1).  Using 
40  bias-estimated  covnriance  lags  from 
200  samples 
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Figure  7.  Spectral  estimate  for  ARMA(3,2).  Using 
40  bias-estimated  covariance  lags  from 
200  samples 


Jf 

f* 

_s 

F.ttci  AMMv) 

ExflcJt.  Wo<U. 

1 

o' 

4 

• 

z 

14 

4  o' 

if  0 

/c 

<L.  KoAlI 

05,4’ 

(5,4)  ) 

i  la 

HC 

HC 

1  Sine,  i.  ^  hloiM, 

0° 

2o 

1 

(5,4’ 

(5,4.) 

4  kll 

He 

uc 

■3  Si«c 

loo 

Zo 

a  h-J.L 

(5,4  ) 

4.  hall 

(2., el 

1 

4  .  Star* 

IUo 

2o 

a 

(2,1  ' 

(  9,8  ) 

4.  Us 

C7o,o  ) 

<.*;<■ ) 

5.  « Mf.(  i,0 

4. 

a.  k.J.L 

C  2,1  1 

(3,.  > 

4  ,  H>U 

(4>,0  ) 

<•«,*  > 

t  A^ma(.s(z) 

Zoo 

4© 

fl  lUai^C 

U,i  ^ 

(4,3) 

k.  **U 

C4^o  ' 

(4,3) 

Table  1.  Experiments.  NC  means  not  computed. 
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